Goldstone modes in Lyapunov spectra of hard sphere systems 



: 

o . 
o ■ 

(N ; 

Oh- 

<: 
o\ : 



Q 

U 
d 



(N 
> 

in 
o 

(N 

m 
o 



> 
X 



Astrid S. de WijrQ and Henk van Beijerer0 

Institute for Theoretical Physics, Utrecht University, 
Leuvenlaan 4, 3584 CE, Utrecht, The Netherlands 
(Dated: February 8, 2008) 

In the study of chaotic behavior, Lyapunov exponents play an important part. In this paper, we 
demonstrate how the Lyapunov exponents close to zero of a system of many hard spheres can be 
described as Goldstone modes, by using a Boltzmann type of approach. At low densities, the correct 
form is found for the wave number dependence of the exponents as well as for the corresponding 
eigenvectors in tangent-space. The predicted values for the Lyapunov exponents belonging to the 
transverse mode are within a few percent of the values found in recent simulations, the propagation 
velocity for the longitudinal mode is within 1 %, but the value for the Lyapunov exponent belonging 
to the longitudinal mode deviates from the simulations by 30%. For higher densities, the predicted 
values deviate more from the values calculated in the simulations. These deviations may be due to 
contributions from ring collisions and similar terms, which, even at low densities, can contribute to 
the leading order. 

PACS numbers: 05.45.Jn 



I. INTRODUCTION 

In the past years, interest in the connections between 
dynamical-systems theory and statistical mechanics has 
increased and many important results have been ob- 
tained. Some of the interest has been directed towards 
the connections between chaoticity and the decay to equi- 
librium. Gallavotti and Cohen 0, put forward a 
chaotic hypothesis, conjecturing that many-particle sys- 
tems as studied by statistical mechanics will generically 
be strongly chaotic. A central role in the study of these 
and related properties is played by the Lyapunov expo- 
nents, which describe the exponential separation or con- 
vergence of nearby trajectories in phase space. 

Many calculations of chaotic properties have been done 
for variations of the Lorentz gas, which is a system of one 
particle bouncing between fixed spherical hard scatterers 
(see, for example, Refs. 0,0,IE@)- Recently, calcula- 
tions have also been done for the largest exponents of 
systems of many freely moving hard spheres 0, Q , which 
is a more realistic model than the single-particle system 
of a Lorentz gas. Simulations for the entire spectrum 
of this system have been done by Posch, Hirschl, and 
Dellago [3, E| ■ The smallest positive and corresponding 
negative exponents from these simulations have received 
a lot of attention because of their unexpected behavior. 
For large enough systems, they are inversely proportional 
to the system length. The tangent-space eigenvectors as- 
sociated with these exponents have a wave like form. 

Attempts have been made to approach these exponents 
by using random- matrix theory by Eckmann and Gat 
[lT| , and by Taniguchi, Dettmann, and Morriss [l2l fl3| . 
Another approach based on kinetic theory has been taken 



by McNamara and Mareschal [Tj] . 

In this paper, we explain some of the behavior of the 
small exponents both qualitatively and quantitatively. 
Sections ITT1 and ITTT1 are a short introduction to Lyapunov 
exponents and the results of the simulations done by 
Posch and Hirschl for hard spheres in two dimensions 
0. In Sec. IIV1 we show that the small exponents are 
in fact due to Goldstone modes. After an explanation of 
the dynamics of hard spheres in Sec. we derive a set 
of equations for the values of the exponents in Sec. IVII 
The equations are derived by using a Boltzmann type 
of approach including a Stofizahlansatz. Finally, we dis- 
cuss the general form of the solutions in Sec. lVlll and the 
quantitative results derived from them in Sec. IVIIII 



II. LYAPUNOV EXPONENTS 

Consider a system with an A/"-dimensional phase space 
r. At time t — the system is assumed to be at an 
initial point 7 in this phase space, evolving with time 
according to 7(70, t). If the system is perturbed by an 
infinitesimal difference S"fo in initial conditions, it evolves 
along an infinitesimally different path 7 + ^7, where ($7 is 
in the tangent space ST. The evolution in tangent space 
is described by 

7(7o + <^7o,0) = 70 + <57o, (1) 
7(7o + ^7o,0 = 7(70,*) +£7(70,*) , (2) 
£7(70,*) = M 7D (i) -<J 7o , (3) 

where M 7o (i) is an A^-dimcnsional matrix, defined by 

^7(70,0 



'70 



(t) 



d-yo 



(4) 
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The Lyapunov exponents are the average rates of growth 
of such infinitesimal changes that are eigenvectors of M 7o , 
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where /J.i(t) is the ith eigenvalue of M l0 (t). In sys- 
tems which are ergodic, almost every trajectory comes 
infinitcsimally close to any point in phase space. This 
means that the Lyapunov exponents are almost indepen- 
dent of the initial conditions. Often the Lyapunov ex- 
ponents are defined not by using M 7o (i), but [M 7o (t) ■ 
M 70 (tY] 1 / 2 . In the latter definition the exponents are 
real. The imaginary components of the Lyapunov expo- 
nents, as we define them here, are also referred to as the 
winding numbers. 

For a classical system of hard spheres without internal 
degrees of freedom, the phase space and tangent space 
may be represented by the positions and velocities of all 
particles and their infinitesimal deviations, 



7, = [r^VjJ , 



(6) 
(7) 



where i runs over all particles and 5"/i is the contribution 
of particle i to Sj. 

In the case of a purely Hamiltonian system, such as 
hard spheres with only the hard particle interaction, the 
dynamics of the system are completely invariant under 
time reversal. Also, due to the incompressibility of flow 
in phase space, the attractor is invariant under time re- 
versal. Therefore, every tangent-space eigenvector that 
grows exponentially in forward time decreases exponen- 
tially in backward time. Since the Lyapunov spectrum 
does not change under time reversal, for every positive 
Lyapunov exponent in such a system there is a negative 
exponent of equal absolute value. This is called the con- 
jugate pairing rule. In systems which are reversible, but 
for which the attractor is not invariant under time re- 
versal, the conditions for and the form of the conjugate 
pairing rule are somewhat different p|. 

Vectors in tangent space which are generated by sym- 
metries of the dynamics of the system do not grow or 
shrink exponentially. They are eigenvectors with Lya- 
punov exponents and are referred to as the zero modes. 
For a system of hard spheres under periodic bound- 
ary conditions, these symmetries and their correspond- 
ing zero modes are uniform translations, Galilei transfor- 
mations, time translations, and velocity scaling. They 
correspond to the initial displacements 
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(Ar o ,0) , 
(0,Av o ) , 
(AtoVi.O) 
(0,AA oVl ) 



= (Ar v ,0) , 
= (0, Av v ) , 



(8) 
(9) 
(10) 

(11) 



where Aro, Avo, Afo, and AAo are constant vectors and 
scalars which are independent of the particle. The quan- 
tities Aro and Avo can have components in all d direc- 
tions of the space. In the case of Galilei transformations 
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FIG. 1: The spectrum of Lyapunov exponents from the 
simulations |g, 1 10L Il5j of 750 hard spheres in two dimen- 
sions at density n = 0.1 in a rectangular box of dimensions 
10 x 75 a' 2 /n, with periodic boundary conditions. Only the 
positive exponents are plotted, since, by the conjugate pair- 
ing rule, the negative spectrum is exactly the opposite. The 
inset shows an enlargement of the bottom right corner. 



and velocity scaling the tangent-space vectors grow lin- 
early, rather than exponentially, with time. Hence the 
corresponding Lyapunov exponents are zero. 



III. LYAPUNOV SPECTRUM OF HARD 
SPHERES 

In principle, M 7o (t) occurring in Eq. © can be cal- 
culated numerically for finite times for any finite sys- 
tem and the eigenvalues can be determined. Posch and 
Hirschl have done molecular dynamics simulations to 
determine the entire Lyapunov spectrum of systems con- 
sisting of many hard disks in rectangular boxes with pe- 
riodic boundary conditions. A spectrum as calculated by 
Posch and Hirschl is displayed in Fig. ^ The eigenvec- 
tors in tangent space belonging to the large exponents 
are typically very localized; only a few particles closely 
together contribute significantly to a given eigenvector. 

When the system is large enough compared to the 
mean free path, a step structure appears in the Lyapunov 
exponents near zero. The size of the steps is inversely 
proportional to the largest dimension of the box. The 
tangent-space eigenvector is distributed over all parti- 
cles, much in the same way as with the zero modes. An 
example of this is shown in Fig. [21 

The tangent-space vectors belonging to the six expo- 
nents in each step appear, on average and to first approx- 
imation, to be linear combinations of the zero modes with 
a sinusoidal modulation. This is very apparent in the ex- 
ample in Fig. [21 The slow modes belonging to a certain 
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FIG. 2: The component of 8v in the short direction 8x is plot- 
ted against the position of the particle in the long direction y 
for A 1493, one of the transverse modes in the first step. These 
data are from the same simulation as the data in Fig. Q The 
corresponding exponent is indicated there with a full box. 

wave vector can be separated into two groups, one con- 
sisting of four longitudinal modes and the other one of 
two transverse modes. The transverse modes are found 
to be linear combinations of sinusoidal modulations of 
the zero modes resulting from a translation or a Galilei 
transformation in the direction perpendicular to the wave 
vector. The longitudinal modes are linear combinations 
of modulations of the four remaining zero modes. The 
transverse modes are non propagating, but the longitudi- 
nal modes propagate through the system. This behavior 
has been confirmed in direct simulation Monte Carlo sim- 
ulations, performed by Forster and Posch Jl5|. For more 
details on these modes, see Refs. [9HIo|. 



FIG. 3: Two particles at a collision in relative phase space. 
The collision normal a is the unit vector pointing from the 
center of one particle to the center of the other. 



time. These modes were first introduced by Goldstonc 
16]. Hydrodynamic modes and phonons in crystals are 
well-known examples. 

In order to calculate the Lyapunov exponents belong- 
ing to these Goldstone modes, one first needs to consider 
the dynamics of the system. 



V. DYNAMICS OF HARD SPHERES IN PHASE 
SPACE AND TANGENT SPACE 

Consider a gas of identical hard spheres or disks of 
diameter a and mass m in d dimensions. The evolution 
in phase space consists of an alternating sequence of free 
flights and collisions. During free flights the particles 
do not interact and the positions grow linearly with the 
velocities, 



IV. GOLDSTONE MODES 

The sinusoidal modes found in the simulations may be 
explained as Goldstone modes. These occur in systems 
with a continuous symmetry, such as the symmetries as- 
sociated with the zero modes. Translation invariance, 
e.g., causes the evolution operator to commute with the 
translation operator, so that they have a set of common 
eigenfunctions. These have the general form 

hi = /k(vj, r y ) exp (ik ■ r;) , (12) 

where the eigenvalues of the operator translating over the 
vector a are of the form exp ik • a. The Goldstonc modes 
are those eigenmodes that for k — > reduce to linear 
combinations of the zero modes. For nonzero values of 
k, they contain a sinusoidal modulation in space of the 
continuous symmetry which grows or shrinks slowly with 



Ti(t) = r I (i ) + (t-<o)v l (t ) , (13) 
v,-(t) = Vi(t ) • (14) 

At a collision, momentum is exchanged between the col- 
liding particles along the collision normal, a = (r,-— Tj)/a, 
as shown in Fig. [3] The other particles do not interact. 
Using primes to denote the coordinates in phase space 
after the collision, we find 

A = r, , (15) 
v- = v< - &(a ■ Vij) , (16) 

where Vjj = Vj — Vj . 

From Eqs. Q and H3I16|1 the dynamics in tangent 
space can be derived Q- During free flight there is no 
interaction between the particles and the components of 
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the tangent-space vector transform according to 

(#)-*(<-«•>■(&). <"» 

/ I (t-t )\ \ 
Z(t-t ) = I | J , (18) 

in which I is the d x d identity matrix. 



At a collision between particles i and j, only the con- 
tributions to the tangent-space vectors of the colliding 
particles are changed '.8|. As shown in Fig. OH an infinites- 
imal difference in the positions of the particles leads to an 
infinitesimal change in the collision normal and collision 
time. The v + <5v are exchanged along a + Sa according 
to Eq. (|16|l . This leads to infinitesimal changes in both 
positions and velocities right after the collision, 
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in which S and Q are the d x d matrices 



A. Equations in fi space 



s = 



(20) 



q _ [(<7 ■ Vij) I + jv.j] • [(» ■ Vjj) I - Vjjoj , 21 . 
a(a-Vij) 

Here the notation ab denotes the standard tensor prod- 
uct of vectors a and b. Let Z(i) be the Nd x Nd matrix 
which performs the transformations of 2(t) on all par- 
ticles. Let L p be the Nd x Nd matrix which performs 
the transformations of Eq. I|19|) on the two particles in- 
volved in collision p and leaves the rest of the particles 
untouched. M 7o (t) in Eq. 10 is a product of these ma- 
trices for the sequence of collisions (1,2, ... ,p) between 
time t and to- Its specific form reads 



'~>o 



(t) = Z(t - t p ) • L p • Z(t p - tp_i)L„_i 



U-Zfa-to) 



(22) 



VI. BOLTZMANN AND ENSKOG EQUATION 

Except in a calculation where the path of every par- 
ticle would be calculated rigorously from the initial con- 
ditions, as done in the molecular dynamics simulations, 
it is impossible to know the matrix M 7o in Eq. J3J ex- 
actly. The tangent-space eigenvector belonging to a given 
Lyapunov exponent will in general depend on the ini- 
tial conditions of all the particles in a much too compli- 
cated way to specify exactly. It is therefore impossible to 
know the tangent-space vector which belongs to any Lya- 
punov exponent exactly. To find the exponents, one has 
to make some statistical approximation without allowing 
contributions along faster growing tangent-space vectors 
to blow up. To this end, we start with assumptions sim- 
ilar to the Stofizahlansatz in the Boltzmann equation. 



To illustrate our calculations we first briefly review the 
Boltzmann and Enskog equations describing the dynam- 
ics of hard-sphere and hard-disk systems at low, respec- 
tively moderate densities. For either one may start from 
the first Bugoliubov-Born-Green-Kirkwood-Yvon hierar- 
chy equation 



df(r,v,t) 
dt 



- v • V r /(r, v, t) + V v • a(r, v, t)f(r, v, t) 



du da na d 1 | a ■ (v — u) | 
/ (2) (r,v',r + aa,u',t) 
-/ (2) (r,v,r - aa,u,t) , 



(T-(V-U)<0 

x 



(23) 



in which denotes the two-particle distribution func- 
tion, u and v respectively u' and v' are the velocities 
before the collision with collision normal a in the direct 
and restituting collisions. 

The second and third term on the left-hand side of the 
equation, respectively, describe the effects of free flight in 
position space and those of the action of external forces. 
The vector a(r, v, t) is the acceleration of a particle due 
to external forces as a function of position, velocity, and 
time. 

In the low-density approximation, Boltzmann's 
StoBzahlansatz approximates the precollisional pair dis- 
tribution functions in this equation by products of one- 
particle distributions. In addition both of these are evalu- 
ated at the same position r, which is allowed if the radius 
a is small compared to the mean free path. The Enskog 
equation is a heuristic generalization of this, known to 
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give a good approximate description of the dynamics up 
to moderate densities (about a quarter of close-packing). 
In this equation the pair-distribution is approximated by 
the product of two one-particle distribution functions, 
evaluated at the actual positions of the two particles, 
and a factor xe, equal to the equilibrium pair correlation 
function at contact between the two particles evaluated 
as a function of the density n((rj + T2)/2) at the point 
halfways between ri and r2 . Notice that this approxima- 
tion becomes exact for a system in homogeneous equi- 
librium. The explicit form of the Enskog equation thus 
becomes 



df(r,v,t) 



dt 



v V r /(r,v,f) 

du da XE{n)na 1 \ a ■ (v — u)| 



5--(v-u)<0 

x [f(r,v',t)f{r + aa,u',t) 
-f(r,v,t)f(r-aa,u,t)] 



(24) 



This equation effectively reduces to the Boltzmann equa- 
tion in the limit n — > 0, when the difference in position 
between the two colliding particles, r.y = aa, may be 
ignored and \e approaches unity. 

In equilibrium the time derivative term vanishes and, 
in absence of external fields, the particles are distributed 
homogeneously This yields the Maxwellian solution 



/(r,v,t) = ncf) M (v) 



d/2 



exp 



(25) 



k B T J * V 2fc B T 
where T is the temperature, related to the average ki- 



netic energy per particle E through E = dk B T '/2. In this 
paper, only the equilibrium system is studied. 

More details on the Boltzmann equation and Enskog's 
theory of dense gases may be found in Refs. 01 an d E3- 



B. Equations in tangent /i space 

To describe the dynamics in tangent space, a general- 
ized Boltzmann equation must be derived for the single- 
particle distribution function in both /i space and "tan- 
gent jU space", /(r,v, Sr, Sv,t). On integration over the 
variables in tangent space, the equation and the solutions 
we are interested in must reduce to Eqs. I|24|) and l|25|l 
respectively. 

For given initial conditions the eigenvectors of M 7o in 
general depend sensitively on the precise values of the 
collision parameters of all collisions, as generated by the 
positions and velocities of all particles. The zero modes 
are exceptions to this. For small k it is to be expected 
that the Goldstone modes behave in a similar way and 
are approximately independent of the collision parame- 
ters of the various collisions. Under those circumstances 
one may expect that the tangent-space vectors Sr and <5v 
can be described by a single-particle distribution function 
that depends smoothly on velocity, position and time, 
just like the velocity distribution in ordinary [i space. If 
in addition one makes the assumption that the distribu- 
tion function of the tangent space vectors of two particles 
about to collide, factorizes in a similar way as the distri- 
bution of their velocities, one ends up with a generalized 
Enskog equation in tangent \i space, which, in absence of 
an external field, is of the form 



J 



df(T,V,dT,6v,t) 



dt 



v • V r /(r, v, Sr, Sv, t) + Sv ■ Va r /(r, v, 5r, Sv, t) 
duda dSs dSuxE (n)na d ~ 1 \ a ■ (v — u)| 



CT-(V-U)<0 

x [/(r, v', Sr', SV, t)f(r + aa, u', 6s', Su' , t) - f(r, v, Sr, 6v, t)f{r - aa, u, Ss, Su, t)] . (26) 



If the tangent /i-space variables Sr and Sv are integrated 
over, this equation reduces to Eq. (|24|l . 

Because <5r and Sv are infinitesimal, the dynamics in 
tangent space are linear in these quantities. Therefore, 
from Eqs. I|19l) and (|26|) one may obtain closed linear 
equations for the time evolution of the average first mo- 
ments (Sr) and (Sv). To this end, multiply both sides of 
Eq. I|26|) by the tangent-space vectors and then integrate 
over them. 



The result is a set of equations for the averages, 

-^<5r(r,v,i) = -v ■ -^-Sr(r,v,t) 

+ 5v(r,v,t) + C s Sr{r,v,t) , (27) 

^-Sv(r,v,t) = v ■ -^-Sv(r,v, t) 
dt or 

+ C s Sv{r,v,t) + C Q Sr(r,v,t) .(28) 
The functions <5r(r, v, t) and 6v(r, v, t) are the averages of 
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the tangent //-space vectors of a particle, as a function of operators Cs and Cq are associated with the matrices S 
its position and velocity, and of time. The linear collision and Q, and given by 



C s 5q(r,v,t) 



C Q <5r(r,v, t) = 



duda XE{n)na \ a ■ (v — u) | (u) 



CT-(V-U)<0 

x { <Sq(r, v', t) + S • [ <Sq(r + a&, u', t) - <5q(r, v', t)] - oq(r, v,t)} , 
d\idaxv(n)na d ~ 1 \ a ■ (v — u)|(/>m(u) 

<T-(V-U)<0 

xQ- [5r(r + aa,u',t) - Sr(r,v',t)] , 



(29) 



(30) 



where <5q can be either or or <5v. In Eq. (|29|l the first two 
terms between braces are gain terms. The last term is 
the loss term. Note that, from Eq. (|21|l . Q is a function of 
the collision parameter and the velocities of the particles 
before the collision. This means that in Eq. l|3Ufl it is 
a function of a, u', and v'. The collision operators are 
proportional to the collision frequency v, which for dilute 
systems is proportional to the number density n. 

C. Fourier transform 

As the translation operators commute with the colli- 
sion operators l|29|l and J30J), solutions to Eqs. (|27|) and 
(|28|l may be found of the form 

<5q(r, v, t) = Aq(v) exp(ik • r + \t) , (31) 

where q is either r or v, and kj = 2nrij/Lj is the jth 
component of the wave vector of the sinusoidal modula- 



tion and A is the Lyapunov exponent. Among these the 
Goldstone modes are those solutions that in the limit of 
vanishing wave number reduce to linear combinations of 
the zero modes. For these modes to stand out among 
the continuum of other modes their wavelength has to be 
large compared to the typical length scale of the mean 
free path, or 

fc|v| < v . (32) 

On substituting Eq. (JSTJ into Eqs. ^ and they 
become eigenvalue equations for the Goldstone modes, 

AAr(v) = -i(k • v)Ar + Av + B$Ar , (33) 
AAv(v) = -i(k-v)Av + B s Av + B Q Ar . (34) 

Spatial propagation, as seen in the simulations for the 
longitudinal modes, may be accounted for by allowing A 
to have an imaginary component. Bs and Bq are the 
Fourier transforms of Cs and Cq, 



B S Aq(v) = 



BqAi-(v) = 



duda \E(n)na \a ■ (v — u)|</>m(u) 



(T-(V-U)<0 

x {Aq(v') + S • [Aq(u') exp (-iak ■ a) - Aq(v')] - Aq(v)} , 
duda-XE(«)na d ~ 1 |(T • (v — u)|</>m(u) 

(T-(V-U)<0 

xQ ■ [Ar(u') exp (-iak • <x) - Ar(v')] , 



(35) 



(36) 



where Aq can be either Ar or Av. One of these, for Eq. i|34|) . This yields 
instance Av, can be eliminated from the equations. One 

must solve Eq. (|3"3f for Av and substitute the result into [(A + ik ■ v — Bs) 2 — Bq]Ar = . (37) 



This equation can be solved by the use of a perturbation 
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expansion in powers of k, provided the mean free path 
is much smaller than the wavelength, as expressed by 
Eq. (32) ■ This is done in Sec. I VII Al 

VII. SOLUTIONS 
A. Perturbation theory 

When Eq. ((32(1 is substituted, one may expand the 
operators and solutions as 

B s = B s 0) +fcB s 1) +fc 2 B^ 2) + ... , (38) 
B Q = B Q 0) +fcBQ 1) +fc 2 B Q 2) + ... , (39) 
Ar = Ar<°> + fcAr^ + fc 2 Ar^ + . . . , (40) 

Note that for k — > 0, the linear operators Bs and Bq 
become identical to Cs and Cq. When acting on linear 
combinations of zero modes Ar' ) they satisfy the prop- 
erties 

Bf Ar(°) = B^ 0) Ar(°) 

= (Ar(°) | B< 0) = (Ar(°) | B Q 0) = , (41) 

where (.|.) represents the inner product defined by in- 
tegration against a Maxwell distribution of the velocity, 
which is the equilibrium distribution. As it turns out, 
B Q 0) has some nontrivial right eigcnfunctions with zero 
eigenvalues, which will have an important effect on the 
limiting values of Lyapunov exponents in the limit of van- 
ishing density. An example of such an eigenfunction is 
Ar(v) = v±h + wiik_L, where v± and v» are the compo- 
nents of v perpendicular and parallel to k. 
In zeroth order Eq. 1)37(1 reduces to 

[(A<°> - Bf ) 2 - B Q 0) ]Ar (0) - . (42) 

The relevant solutions to this are the zero modes, with 
A<°> = 0. This means that the Goldstone modes to lead- 
ing order in k are the zero modes with a sinusoidal modu- 
lation, in nice agreement with the findings in Refs. [9l[To|. 
In linear order, one finds with the aid of Eq. I|41|) 

[-Bf (ik-v-B^J-B^lArPO 

= -[(Bf) 2 -B^]Ar«, (43) 

where k is the unit vector in the direction of k. This is 
an equation for Ar' 1 '. Its formal solution is 

Ar« = [(B< 0) ) 2 - BgV 

[B^ 0) (tk • v - B^) + B Q 1} ] Ar^ . (44) 

This form suggests that Ar' 1 ' is of zeroth order in n, just 
as Ar' ), but this is actually not the case, because the 
operator [(Bg°' ) ) 2 — Bq acts on functions with non van- 
ishing components along the nontrivial right zero eigen- 
functions of Bq"*. This yields contributions to Ar' 1 ) of 



order 1/n. One might wonder whether this could cause 
a divergence in the limit of vanishing density, but that is 
not the case because of the restriction imposed on k by 
Eq. 

The second order equation involves the first-order Lya- 
punov exponent A^, the second-order Lyapunov ex- 
ponent A^ 2 ' , and the second-order tangent-space vector 

Ar (2) 7 

+ (A( 1 )+zk-v-B( 1) ) 2 -B Q 2) ]Ar(°> 
+ [{-Bf,zk.v-BW} + -B«]Ar« 
+ [(Bf) 2 -B Q 0) ]Ar( 2 )-0, (45) 

where {., .}+ is used to denote the anticommutator of two 
operators. On taking the inner product with Ar' ), all 
terms involving Ar^ 2 ) vanish as a consequence of Eq. (f 4 1 1> . 
The resulting set of equations reads 

(Ar (0) |[(A (1) + zk • v - B^) 2 - B Q 2) ]|Ar(°)) 

+ (Ar<°>|[-(ik • v - B^)Bf - B^IArW) 

= . (46) 

Since Ar*- ) is a linear combination of three independent 
zero modes, Eq. ((46(1 actually has to be read as a 3 x 3 ma- 
trix equation involving the matrix elements between the 
various zero modes. In principle all of these are second- 
order polynomials in A*- 1 ). The eigenvalues, as usual, fol- 
low from the condition that the determinant of the matrix 
vanishes as a function of A^ . 



B. General form of the solutions 

In order to investigate the general structure of Eq. 1(46(1 
it is useful to organize the zero modes for Ar^ ^ as 

Arf=k ± ; Arf^k; Ar[»)=^fv, (47) 

where (3 = l/(k-sT). The first mode consists of a perpen- 
dicular displacement, i.e., a spatial translation normal to 
the wave vector, the second mode to a parallel displace- 
ment, and the third one to a time translation. 

The first mode is odd in kj^ and the last two even; the 
first two modes are even in v and the last one odd. The 
collision operators B$ and Bq as well as the function k • v 
are odd in kj^ to every order. The operators Bg™" 1 and Bq 1 " 1 
are even in v for even n and odd for odd n. On the basis 
of these parity properties it follows immediately that the 
structure of Eq. 146(1 , written as a matrix equation on the 
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basis (|47|l . is restricted to 





(1 









(A«) 2 





1 


i 

















(48) 



The constants x and y are determined by temperature 
and by the form of the collision operators. From this 
it becomes clear that the equation can be split into two 
parts, one for the perpendicular zero mode, the trans- 
verse part, and one for the parallel zero mode and the 
time mode, the longitudinal part. From the general form 
of the matrices one can derive the general form of the 
Lyapunov exponents A = fcA^ 1 -* to be 



Atrans 
Along 



(49) 
(50) 



where y\ and yi are functions of x V) |i , xi| )V , j/iin, and y v ,v 
If y±,± > 0, the Lyapunov exponent of the transverse 
mode is real and therefore the mode is of the same form as 
in the simulations reported in 0,^3- If 2/2 > 0, the lon- 
gitudinal Lyapunov exponents have both real and imagi- 
nary components, and these modes also have the form of 
the longitudinal modes found in the simulations. 



C. Density expansion 

The Stofizahlansatz is an approximation that for many 
purposes, e.g., the derivation of hydrodynamic equations 
from the Boltzmann equation with explicit expressions 
for the transport coefficients 0], becomes exact in the 
limit of vanishing density. Therefore we want to inves- 
tigate the behavior of our equations in this limit and 
compare to the results found in the simulations. In the 
limit of vanishing density Eq. (|46|l becomes 

(ArW|[A (1) +z(k- v)] 2 |ArW) 

- (Ar (0) |[z(k • v)B< 0) + B^]|Ar (1) ) = . (51) 

For this equation it is crucial indeed that Ar^ is of the 



order of n 
tions of B 



. If there were no nontrivial right eigenfunc- 
with eigenvalue 0, Ar^ would be one order 
of n higher, and the second term would not contribute to 
the Lyapunov exponents in the limit of vanishing density. 

In the following section we will further discuss the ac- 
tual magnitudes of the two terms in Eq. I|51(l . 



VIII. RESULTS AND DISCUSSION 

If only the first term in Eq. I|46|l is kept, the calculation 
is fairly simple. From now on we choose d = 2. The same 



calculations can easily be repeated for d = 3, but there 
are far fewer simulation results to compare to. Equation 
(I48JI in this approximation becomes 



( A (l))2/^ 



/ 1 



\0 1 




independent of the density. The solutions for the Lya- 
punov exponents then are 



= ±- 



k 



» (±0.978 ±i 0.676)^ 



(53) 
(54) 
(55) 



The structure of the corresponding eigenvectors is indeed 
like that found in simulations 0, flol | . 

To calculate the contribution from the second term in 
Eq. i|46|) to the leading order of the Lyapunov exponents, 
one has to choose a suitable basis in which to express 
the function ArW(v). The basis must be orthogonal 
with regard to the chosen inner product (.|.). Next, the 

matrix elements of the operators Bg and 
calculated between elements of the basis. 

A simple, but suitable, basis is the set of functions 
that are products of Hcrmitc polynomials in the compo- 
nents of v-y/m/3/2 parallel and perpendicular to the wave 
vector k. The Hermite polynomials Hi(x) form a com- 
plete orthogonal basis with regard to integration against 
exp (— x 2 ), and therefore their products will be orthogo- 
nal under the inner product used here. The solution to 
Eq. I|37|) can thus be expanded as 



Bq^ must be 



Ar ( V ) = ^2 e l c Lp,q H p( v \\) H q( v ±) . 



(56) 



where I can be either _L or 



X is kj_ and e || is k. 

By truncating all expressions at some finite order in 
the polynomial expansion, one finds approximate values 
for A^ 1 ). For good convergence one has to go beyond the 
zeroth- and first-order Hermite polynomials. In the Ap- 
pendix more details are given on the matrix representa- 
tions of the truncated operators and on the convergence 
of the Lyapunov exponents in dependence on the order 
of truncation. 

To sixth order in the polynomial expansion, the results 
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FIG. 4: The Lyapunov exponents for transverse and longi- 
tudinal modes from the simulations 0, Ho|l compared to the 
present calculations. 
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FIG. 5: The velocities of the longitudinal mode from simu- 
lations compared to the present calculations. 

for n — > and d = 2 are 



Atrans — 
'Mong — 
^long — 




(57) 
(58) 
(59) 



With respect to Eqs. I|53|) - (|55|l the corrections are largest 
for the Lyapunov exponent of the longitudinal mode. 

For low densities the form of the modes is predicted 
correctly by the calculations; the modes are split into 



non propagating transverse and propagating longitudi- 
nal modes. For number density n = pa 2 = 0.02, the 
Lyapunov exponents from the simulations are 



-^trans 
^long 
^long 




(60) 
(61) 
(62) 



The calculated Lyapunov exponent of the transverse 
mode and the propagation speed of the longitudinal mode 
compare to the values from the simulations, within 2% at 
n = 0.02. The Lyapunov exponent for the longitudinal 
mode deviates by about 30%. 

The results for higher densities are displayed in Fig.0] 
With increasing density the calculated values deviate in- 
creasingly from the simulation results. For the longitudi- 
nal mode the predicted real part of the Lyapunov expo- 
nent even drops to and the exponent becomes purely 
imaginary. 

The deviations from the simulations can be attributed 
to contributions from ring terms and possibly other con- 
tributions to a generalized Bq operator that are at most 
of order n 2 . From Eq. I|44|l one sees that such terms, 

working on the non trivial zero eigenfunctions of Bq\ 
contribute to the leading order terms in the density ex- 
pansion of Ar (1) , just like (B^ 0) ) 2 . Therefore they have to 
be included in the second term in Eq. (51). So in contrast 
to usual applications of kinetic theory, where ring terms 
only contribute to higher orders in the density, in the 
present case they contribute to the leading order. In the 
present calculation these contributions are not included, 
but we are actively working on their evaluation. Simi- 
larly, the ring terms will contribute to higher orders in 
a density expansion of the Lyapunov exponents. These 
contributions may be responsible for the discrepancies 
between simulation results and Enskog theory for higher 
densities, which show up in Figs. E| and For more de- 
tails on ring terms in kinetic theory, see Ref. [l9| . 

It is interesting to compare our results to those by Mc- 
Namara and Mareschal [lj| , who also based their work on 
kinetic theory calculations. They do not derive equations 
for the distribution functions, but go directly to hydro- 
dynamic like equations for the moments. To close these, 
they make hypotheses to factorize the fluxes. The result- 
ing values for the Lyapunov exponents in the low-density 
limit are less close to the simulation values than those 
from our calculations. It is not clear that in this treat- 
ment the effects of the non trivial zero eigenfunctions of 
Bq-' are accounted for. 

Forster and Posch have also done simulations on sim- 
ilar systems with soft potentials |2jj. They roughly find 
a branch again of Lyapunov exponents close to zero, but 
the sinuoidal structure of the corresponding modes is 
much less clear 20] . It would be very interesting to cal- 
culate the Lyapunov exponents with kinetic theory meth- 
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ods also for this case. It would also be very interesting 
to look at small Lyapunov exponents in non equilibrium 
systems. However, in such systems the calculations be- 
come more complicated because the stationary velocity 
distributions are not Maxwellian any more. 



IX. CONCLUSION 

In this paper, we have demonstrated how Lyapunov ex- 
ponents close to zero can be related to Goldstone modes. 
We found the correct types of behavior in dependence 
on the wave number of the exponents and their tangent- 
space eigenvectors. This was achieved through a kinetic 
theory approach, in which we used a molecular chaos as- 
sumption for the pair distribution function to derive an 
equation similar to the Enskog equation. For low densi- 
ties this reduces effectively to a generalized Boltzmann 
equation. 

The calculated values for the exponents belonging to 
the transverse modes at low densities are within a few 
percent of the values found in the simulations |9|,[nj- The 
propagation velocity for the transverse mode is within 1% 
of the simulation values. The value for the Lyapunov ex- 
ponent belonging to the longitudinal mode deviates from 
the simulations by 30%. For higher densities, the pre- 
dicted values deviate increasingly more from the values 
found in the simulations. These deviations are proba- 
bly due to contributions from ring collisions and similar 
terms. In most applications of the Boltzmann equation 
and the StoBzahlansatz such terms produce contributions 
to the relevant quantities which are one order of higher 
order in the density, but in the problem at hand they 
turn out to contribute to the leading order. 



the zero modes are 



(0) 



(1,0,0,0,0,0), (Al) 
(0,1,0,0,0,0), (A2) 
(0,0 ±\/2,0,0 ,|V2) . (A3) 



Here, the subscript 1 indicates that basis functions up to 
first order in v have been included. From Eqs. (|43|l and 
(|46|l it follows that the operator B$ described in Eq. I|35|) 
is only needed up to first order in k. One finds for the 
matrix elements of B$ in the expansion of Eq. (|38|) 
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(A4) 



(A5) 



The first three contributions to Bq as expanded in 
Eqs. H36(l and (|39J) have similar matrix representations 
of the form 
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APPENDIX A: EXPANSION IN HERMITE 
POLYNOMIALS 

In order to solve Eq. I|46l) , one must write the operators 
in Eqs. (|35|l and (|36|l as matrices between the basis func- 
tions described in Eq. (J5BJ. From now on we take d = 2, 
but the same calculations can easily be done for three di- 
mensions. We only show results for basis functions of up 
to linear order in v. In Eq. I)56|l p and q can be equal to 
or 1 . In fact one has to include higher powers to find good 
approximations for the solutions to the original equa- 
tions. If the first component is the component parallel 
to k, the basis is ordered as (1,0); (0,1); (y / (3m/2v^0); 

(0,y^m^«||); (y/0m/2v ± ,O); (0,^/^/2^). All co- 
efficients are given to leading order in n. In this notation 
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(A6) 



(A7) 



(A8) 
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The operators ik • v and — (k • v) 2 can also be written in 
this way. One finds 




-(k- 



/0 i 0\ 

i 

1 
z 


Voooooo/ 

/ 1 \ 

1 

3 

3 

1 

V 1 / 



(A9) 



(A10) 



With these matrices and Eq. (|43|l the vectors for Ar^ 
may be expressed in terms of Ar^ ^ up to first order in the 
polynomial expansion in v. With the orthogonality rela- 
tion between Ar^ and Ar^ ^ , mentioned in Sec. IVII Al 
this yields 



Ar« = 



2i 



9nay / 7r 



fo \ 



1 
1 
1 

V-l / 



Ar 



(0) 



(All) 



These matrices can be used to find the 3x3 matrices in 
Eqs. H46JI and l)48[l. The equation to leading order in n 
then becomes 




(A12) 



Here, the indices of the matrix on the left hand side 
are ordered according to , Ar ^ , Ari°' ) . The matrix 
can be factorized into two parts. One part describes the 
transverse mode and produces a simple quadratic equa- 
tion for A^ 1 - 1 , with the solution 



A ( 1} ^pm = ± i ^ ±Q 62361 _ (A13) 

The remaining part yields the longitudinal mode. The 
solutions to the fourth order equation for A^ 1 ) are 

ps ±0.639552 ±0.493231 i . (A14) 

The same calculation can be done with larger subsets of 
the basis. The results are shown in Table [I] 

Using functions up to an odd power in v is different 
from using functions up to an even power, because the 
odd powered functions contribute to different matrix el- 
ements than the even powered functions. To determine 
whether the solutions have converged, one must therefore 
look at the behavior as the maximum power is increased 



by steps of 2. The error in the results using up to sixth 
powers in v can be estimated by comparing the values 
with the results for powers in v up to four. The error 



n 


transverse A^ 1 ' \J 4p 


longitudinal 


1 


±0.62361 


±0.639552 ± 0.463231 i 


2 


±0.62361 


±0.422807 ± 0.499026 i 


3 


±0.626194 


±0.424806 ± 0.499105 i 


4 


±0.626194 


±0.428599 ± 0.498952 i 


5 


±0.626254 


±0.428645 ± 0.498954 i 


6 


±0.626254 


±0.429104 ± 0.498953 i 



TABLE I: The Lyapunov exponents and the propagation 
velocities for the longitudinal mode calculated using products 
of Hermite polynomials in v up to different orders. 



in the solutions when using up to sixth powers of v in 
the basis functions appears not to be much larger than a 
tenth of a percent, except in the case of the longitudinal 
mode, where it might be of the order of a few tenth of a 
percent. 
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